根据16S预测微生物群落功能最全攻略
点击上方蓝色「宏基因组」关注我们！专业干货每日推送！
声明：本文转载自“微生物生态”公众号，点我或文末阅读原文跳转原文。“微生物生态”，干货多，同行值得关注。本文是此类分析工具的简介和选择指南，具体的网页和本地命令行操作，会在将来的文章中详细介绍。
最近，越来越多的证据表明：自然环境（如海洋、土壤等）中的微生物群落功能（functional）组成而非物种（taxonomic）组成与环境因子密切相关；换言之，相似环境中的微生物群落功能更相似，而行使功能的微生物物种组成可能差异较大（Gibbons et al. 2017; Louca et al. 2016; Nelson et al. 2016）。这说明，除了揭示环境中有哪些微生物之外，揭示微生物群落功能轮廓尤为重要。目前，微生物生态研究中常用的揭示微生物群落功能的方法有宏基因组测序、宏转录组测序、宏蛋白组测序、宏代谢组分析等。这些方法优点突出，能较准确、真实地反映不同层面的微生物群落功能特征。但其价格较高，一般实验室难以承受大批量样本的实验；另外，数据量巨大，数据处理也是难点。有没有较经济、又适合大样本的方法呢？有！那就是基于marker基因扩增子高通量测序的功能预测。上次卢瑟菌给大家介绍了根据真菌ITS序列预测真菌群落功能的工具——FUNGuild（点击了解更多）
今天，卢瑟菌就和大家介绍基于原核16SrDNA高通量测序结果对微生物群落功能（function）或表型（phenotype）进行预测的四种方法——PICRUSt、Tax4Fun、FAPROTAX及BugBase。
简介
PICRUSt全称为Phylogenetic Investigationof Communities by Reconstruction of Unobserved States，由Langille等人于2013年开发，文章发表在Nature Biotechnology上（Langille et al. 2013）。它是最早被开发的基于16S rRNA基因序列预测微生物群落功能的工具，包括在线版（http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=PICRUSt_normalize）和基于MacOS X或Linux系统的下载安装版（http://picrust.github.io/picrust/install.html#install）。
原理
如图1所示，其预测过程分两步：
（1）基因内容预测（gene content inference）。该步先对Greengenes数据库的“closed reference”序列划分OTU后构建进化树，通过祖先状态重构（Ancestralstate reconstruction）算法并结合IMG/M数据库，预测出树中未进行全基因组测序OTU的基因组信息。
（2）宏基因组预测（metagenome inference）。将16SrDNA测序结果与Greengenes数据库进行比对，挑选出与“closed reference”数据库相似性高的（默认为≥97%）OTU；根据OTU对应基因组中16SrDNA的拷贝数信息，将每个OTU对应序列数除以其16S拷贝数来进行标准化；最后，将标准化的数据乘以其对应的基因组中基因含量从而实现宏基因组预测的目的。获得的预测结果可以通过KEGG Orthology、COGs或Pfams等对基因家族进行分类。
图1 PICRUSt预测微生物功能原理流程图（Langille et al. 2013）
特点
（1）因开发者当初开发PICRUSt时基于Greengenes数据库为参考，用该分析时只识别16S序列与Greengenes数据库（GG13.5或GG12.5）比对的结果；
（2）预测结果的准确度可通过参数NSTI（Nearest Sequenced Taxon Index）值的大小来评估。NSTI表征某样品中所有微生物OTU与其亲缘关系最近的已测序基因组间系统进化距离的平均值，因此，该值越小表示预测结果越可信。如图2，通过对不同环境来源的微生物群落功能预测结果及其宏基因组数据比较后发现：PICRUSt对于人肠道微生物样品的预测结果最好（平均NSTI =0.03 ± 0.02 s.d.），其次是土壤样品（平均NSTI= 0.17 ± 0.02 s.d.），其他哺乳动物肠道样品波动较大（NSTI = 0.14 ± 0.06 s.d.），而对于研究较少的高盐微生物席样品的预测准确度最低（NSTI= 0.23 ± 0.07 s.d.）。
图2 不同环境微生物组PICRUSt预测准确度比较（Langille et al. 2013）
（3）对于低16S测序深度的样品也可获得较准确的预测结果。通过比较PICRUSt预测与宏基因组测序随测序深度变化的结果准确度后发现，即使低16S测序量（105条序列）也可获得准确度较高的功能预测结果（图3）。这说明PICRUSt不仅适用于高通量测序后的功能预测，对于传统克隆文库/Sanger测序获得的低测序深度的样品也有较好的预测准确度。
图3 随测序深度增加PICRUSt预测与宏基因组测序结果准确度比较（Langille et al. 2013）
（4）无论对古菌还是细菌都有较高的预测准确度。如图4，无论细菌（mean= 0.95 ± 0.05 s.d., n = 2,487）还是古菌（mean = 0.94 ± 0.04s.d., n = 103），PICRUSt都有较高的预测准确度。
图4 细菌和古菌基因组进化树中PICRUSt预测准确度情况（Langille et al. 2013）
（5）对不同功能基因分组的预测准确度存在差异，但均在可信范围。如图5，对于核心功能基因分组，如遗传信息加工的基因家族预测准确度最高（0.99± 0.03 s.d.）；而对于那些在不同基因组间易变和易发生基因横向转移的基因家族，如环境信息加工基因家族的预测准确度相对较低（0.95 ± 0.04 s.d.）；但整体而言，预测的准确度都高于0.82，说明预测结果可靠。
图5 PICRUSt对基因组功能模块预测准确度差异（Langille et al. 2013）
示例
基于16S rDNA序列的PICRUSt功能预测，可获得不同level（1~3）的KEGG Pathways的功能基因丰度富集情况（图6），也可获得具体与某个功能过程（如硫代谢等）相关的基因丰度情况（图7）。
图6 通过PICRUSt功能基因预测得到的KEGG通路中的基因相对丰度示例（Sun et al. 2016）
图中展示了西藏猕猴在不同季节（winter和spring）肠道微生物的KEGG功能通路预测结果，*表示在不同季节差异显著的功能富集。
图7 通过PICRUSt预测得到的硫代谢相关基因的丰度示例（Wang et al. 2016）
图中显示了不同地点沉积物原核微生物群落硫代谢相关功能基因的PICRUSt预测结果，右侧红色字体表示异化硫酸盐还原相关的基因，蓝色表示同化硫酸盐还原相关基因。
简介
Tax4Fun是Aßhauer等人2015年开发的通过16S高通量测序数据预测微生物群落功能的方法，文章发表在Bioinformatics上（Aßhauer et al. 2015）。其设计思路与PICRUSt类似。Tax4Fun无网页在线版，只能基于R包进行线下分析。相关资料下载链接：http://tax4fun.gobics.de/。
原理
首先，将16S高通量测序数据通过QIIME或SILVAngs平台，基于SILVA数据库对OTU进行物种分类；基于分类结果，再根据NCBI的基因组注释对16S拷贝数进行标准化；最后通过构建SILVA分类与KEGG数据库中原核分类间的线性关系，实现对微生物群落功能的预测。
特点
据Tax4Fun文章介绍，其最大优点是对微生物群落功能预测的准确性要比PICRUSt好。作者对不同环境来源的样品（人类肠道、哺乳动物肠道、高盐微生物席及土壤样品）分别通过Tax4Fun和PICRUSt进行功能预测，并与相应的宏基因组结果进行相关性分析后发现：基于QIIME或SILVAngs平台结合Tax4Fun的预测结果在所有测试样品中与实际测得宏基因组结果的相关性均好于PICRUSt的预测结果，即Tax4Fun对微生物群落功能的预测准确性要优于PICRUSt（图8）。卢瑟菌个人认为，一方面，从预测原理来讲，PICRUSt中相当比例OTU的基因组是经祖先状态重构算法预测出来的，并非真实的基因组信息，而Tax4Fun都是基于KEGG库中已测序注释的原核基因组信息，这可能是决定两者准确度差异的主要原因；另一方面，Tax4Fun选择的是SILVA数据库，而PICRUSt是Greengenes数据库，前者的更新速度更快。
图8 不同样本宏基因组及Tax4Fun、PICRUSt预测结果相关性比较（Aßhauer et al. 2015）
示例
与PICRUSt类似，通过Tax4Fun可预测出不同level（1~3）的KEGG Pathways功能基因丰度富集情况，也可获得具体与某个功能过程（如硫代谢等）相关的基因丰度状况（示例参见PICRUSt示例）。
简介
FAPROTAX取词自Functional Annotation of Prokaryotic Taxa，是Louca等人为解析微生物群落功能于2016年创建的基于原核微生物分类的功能注释数据库，文章发表在2016年的Science上（Louca et al. 2016）。FAPROTAX是基于目前对可培养菌的文献资料手动整理的原核功能注释数据库，其包含了收集自4600多个原核微生物的80多个功能分组（如硝酸盐呼吸、产甲烷、发酵、植物病原等）的7600多条功能注释信息。作者编写了一套python脚本来运行预测，输入文件格式可以是SILVA或Greengenes数据库生成的OTU分类表或BIOM文件。相关资料下载地址：http://www.zoology.ubc.ca/louca/FAPROTAX/lib/php/index.php?section=Home。
原理
FAPROTAX原理与以前提到的真菌功能预测FUNGuild类似。如图9，作者先根据文献资料（Begrey's Manualof Systematic Bacteriology、The Prokaryotes、The International Journal of Systematic Bacteriology等）手动构建了联系物种分类与功能注释的FAPROTAX数据库；后又编写了联系OTU分类表与FAPROTAX数据库的python脚本；最后，只要将基于16S的OTU分类表通过python脚本就可以输出微生物群落功能注释预测结果。
图9 FAPROTAX预测微生物群落功能原理示意图（Louca et al. 2016）
特点
FAPROTAX较适用于对环境样本（如海洋、湖泊等）的生物地球化学循环过程（特别是碳、氢、氮、磷、硫等元素循环）进行功能注释预测。因其基于已发表验证的可培养菌文献，其预测准确度可能较好，但相比于上述PICRUSt和Tax4Fun来说预测的覆盖度可能会降低。与PICRUSt和Tax4Fun类似，FAPROTAX依赖于16S序列的分类结果，较好的分类结果（能分辨到属种水平的物种比例较高）才能得到较好的预测结果。预测结果中可能出现一个OTU对应多个功能分组的情况，如被注释到硫磺单胞菌属（Sulfurospirillum）的OTU即可进行硝酸盐氨化，又可以进行砷酸呼吸和发酵。
示例
FAPROTAX可根据16S序列的分类注释结果对微生物群落功能（特别是生物地化循环相关）进行注释预测（图10）。
图10 FAPROTAX预测的不同分层海水样品微生物群落功能轮廓（Louca et al. 2016）
图中横轴所示为表层（SRF，surfacewater）、叶绿素含量最高层（DCM，deepchlorophyll maximum）、混合层（MIX，mixedlayer）和海洋中层（MES，mesopelagic）海水样品中微生物群落功能轮廓，纵轴所示为包括碳、氢、氮、硫等元素循环相关及其他诸多功能分组。颜色越深代表样品中该类群微生物的相对丰度越高。
简介
BugBase是Dan Knights课题组开发的用于对微生物组数据进行高水平表型（high-level phenotypes）分类的工具，目前相关文章正在整理之中，但该工具已开放，可免费使用（bugbase.cs.umn.edu/index.html）。该工具可对微生物群落根据七类表型进行分类：革兰氏阳性（Gram Positive）、革兰氏阴性（Gram Negative）、生物膜形成（Biofilm Forming）、致病性（Pathogenic）、移动元件（Mobile Element Containing）、氧需求（Oxygen Utilizing，包括Aerobic、Anaerobic、facultatively anaerobic）及氧化胁迫耐受（Oxidative Stress Tolerant）。输入由Greengenes数据库分类后的OTU表格（BIOM格式），即可快速实现对上述表型的分类预测。若同时输入Mapping文件还可以实现对分组变量的作图及统计比较分析。目前有在线网页版（bugbase.cs.umn.edu/upload.html，数据<15M适用）和线下安装版（bugbase.cs.umn.edu/downloads.html，>15M适用）。
原理
由于文章还未发表，因此具体原理未给出。但推测基本原理与上述三种类似，先根据Greengenes等数据库对16S数据进行物种分类，再根据分类结果结合KEGG数据库中原核功能基因的注释或根据文献对可培养菌的报道实现表型预测。
示例
通过BugBase可快速地对16S高通量结果根据七类表型进行分类比较（图11）。
图11 通过BugBase预测得到的细菌群落表型分析示例（Thomas et al. 2016）
图中展示了结肠癌患者及健康人群的结肠样本细菌群落在好氧、厌氧、革兰氏阴性、阳性、生物膜形成及兼性厌氧等表型上的相对丰度差异，p值经Wilcoxonrank sum test计算得到，p <0.05表示差异显著，三条线自上而下分别表示上四分位、平均值及下四分位。
至此，目前常用的基于微生物marker基因的微生物群落功能预测工具介绍完毕。表1从不同角度列出了这些工具的特点，以方便大家选择使用。另外，根据卢瑟菌的个人理解，简单总结各工具的选择策略，仅供参考：
①若想获得原核微生物的KEGG代谢通路或相关酶等类似宏基因组测序获得的功能预测信息，选择PICRUSt和Tax4Fun。其中前者有定量的预测准确度评估参数NSTI，若NSTI≤0.17则说明PICRUSt的预测准确度较好。虽然Tax4Fun文中指出其预测准确性高于PICRUSt，但缺乏类似于NSTI的定量参数评估。因此，当NSTI >0.17时，建议两种方法都试一下，看看结果是否一致。另外，除了考虑预测的准确度之外，预测力度或覆盖度（即能被预测的OTU的序列数占总序列数的比例）也是要考虑的参数。
②想了解原核微生物群落生态功能，尤其是在碳、氢、氮、磷、硫等元素的生物地球化学循环中的功能预测情况，选择FAPROTAX。
③想获得原核微生物群落的好氧、厌氧、兼性厌氧、生物膜形成、革兰氏阳性/阴性等高水平的表型分类情况，选择BugBase。
④想获得真菌群落的共生、腐生、致病性等生态型分类，选择FUNGuild。
⑤上述所有基于marker基因的功能预测方法，都不能完全代表真实的微生物群落功能，若想严谨地证明微生物群落确实存在某功能（基因），仍需要结合宏基因组、宏转录组、宏蛋白组、宏代谢组及其他微生物功能、活性测定方法来共同佐证。
表1 常用微生物群落功能预测工具比较
参考文献
Aßhauer K P,Wemheuer B, Daniel R, et al. Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data[J]. Bioinformatics, 2015, 31(17): 2882-2884.
Gibbons S M. Microbial community ecology: Function over phylogeny[J]. NatureEcology & Evolution, 2017, 1: 0032.
Langille M G I, Zaneveld J, Caporaso J G, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences[J]. Naturebiotechnology, 2013, 31(9): 814-821.
Louca S,Parfrey L W, Doebeli M. Decoupling function and taxonomy in the global oceanmicrobiome[J]. Science, 2016, 353(6305): 1272-1277.
Nelson M B, Martiny A C, Martiny J B H. Global biogeography of microbial nitrogen- cycling traits in soil[J]. Proceedingsof the National Academy of Sciences, 2016, 113(29): 8033-8040.
Sun B, Wang X, Bernstein S, et al. Marked variation between winter and spring gutmicrobiota in free-ranging Tibetan Macaques (Macaca thibetana) [J]. Scientificreports, 2016, 6.
Thomas A M,Jesus E C, Lopes A, et al. Tissue-associated bacterial alterations in rectalcarcinoma patients revealed by 16S rRNA community profiling[J]. Frontiersin Cellular and Infection Microbiology, 2016, 6.
Wang K, Ye X, Zhang H, et al. Regional variations in the diversity and predictedmetabolic potential of benthic prokaryotes in coastalnorthern Zhejiang, East China Sea[J]. Scientific Reports, 2016, 6: 38709.
更多相关文章，请点击“宏基因组”文章目录阅读；
想了解更多16S扩增子、宏基因组文献和分析相关文章，快关注“宏基因组”公众号，干货每日推送
系统学习生物信息，快关注“生信宝典”